Analysis of DESeq2 results with gene set enrichment

Analysis of DESeq2 results with gene set enrichment

0 Boilerplate

Libraries and external code needed:

Set working directory where the script is

if (Sys.getenv("USER")=="JQ") {
  setwd("/Users/JQ/Documents/_CODE REPOS/GitHub/Da_RNAseq")
} else if (Sys.getenv("RSTUDIO")==1) {
  setwd( dirname(rstudioapi::getSourceEditorContext(id = NULL)$path) ) # gets what is in the editor
} else {
  setwd(here::here())
  d <- str_split(getwd(),'/')[[1]][length(str_split(getwd(),'/')[[1]])]
  if (d != 'Da_RNAseq') { stop(
    paste0("Could not set working directory automatically to where this",
           " script resides.\nPlease do `setwd()` manually"))
    }
}
getwd()

## [1] "/Users/JQ/Documents/_CODE REPOS/GitHub/Da_RNAseq"

Path to definitive images (outside repo):

figdir <- paste0(c(head(str_split(getwd(),'/')[[1]],-1),
                   paste0(tail(str_split(getwd(),'/')[[1]],1), '_figures')),
                 collapse='/')
dir.create(figdir, showWarnings = FALSE)

1 Getting ready

This gets us the DGE data from DESeq2, identified by FlyBase/Ensembl ID and gene symbol:

# experimental design and labels
targets <- readRDS('output/targets.RDS')
# DEG data
DaKD_deg <- readRDS('output/Control_vs_DaKD.RDS')
DaOE_deg <- readRDS('output/Control_vs_DaOE.RDS')
DaDaOE_deg <- readRDS('output/Control_vs_DaDaOE.RDS')
ScOE_deg <- readRDS('output/Control_vs_ScOE.RDS')
head(
  ScOE_deg %>% dplyr::select(gene_symbol, ensemblGeneID, baseMean, log2FoldChange, padj),
  1
)

2 Gene set enrichment analysis of custom gene sets

2.1 Cell-type specific transcriptomes (Dutta et al., 2015)

Run GSEA

First let’s run clusterProfiler::GSEA:

dutta_gmx <- import_from_gmx(paste0(getwd(),'/resources/Dutta.gmx'))

daKD_gse_dutta <- GSEA(geneList=make_degrank(DaKD_deg, mode='log2fc', key='ensemblGeneID'),
                       exponent = 1,
                       minGSSize = 1,
                       maxGSSize = 1000,
                       eps = 0,
                       pvalueCutoff = 1,
                       pAdjustMethod = "BH",
                       TERM2GENE = dutta_gmx,
                       TERM2NAME = NA,
                       verbose = TRUE,
                       seed = FALSE,
                       by = "fgsea")

daOE_gse_dutta <- GSEA(geneList=make_degrank(DaOE_deg, mode='log2fc', key='ensemblGeneID'),
                       exponent = 1,
                       minGSSize = 1,
                       maxGSSize = 1000,
                       eps = 0,
                       pvalueCutoff = 1,
                       pAdjustMethod = "BH",
                       TERM2GENE = dutta_gmx,
                       TERM2NAME = NA,
                       verbose = TRUE,
                       seed = FALSE,
                       by = "fgsea")

dadaOE_gse_dutta <- GSEA(geneList=make_degrank(DaDaOE_deg, mode='log2fc', key='ensemblGeneID'),
                       exponent = 1,
                       minGSSize = 1,
                       maxGSSize = 1000,
                       eps = 0,
                       pvalueCutoff = 1,
                       pAdjustMethod = "BH",
                       TERM2GENE = dutta_gmx,
                       TERM2NAME = NA,
                       verbose = TRUE,
                       seed = FALSE,
                       by = "fgsea")

scOE_gse_dutta <- GSEA(geneList=make_degrank(ScOE_deg, mode='log2fc', key='ensemblGeneID'),
                       exponent = 1,
                       minGSSize = 1,
                       maxGSSize = 1000,
                       eps = 0,
                       pvalueCutoff = 1,
                       pAdjustMethod = "BH",
                       TERM2GENE = dutta_gmx,
                       TERM2NAME = NA,
                       verbose = TRUE,
                       seed = FALSE,
                       by = "fgsea")

GSEA plots

Now I need to make GSEA plots from the results. I will use genekitr for this.

# the gsea objects generated from clusterProfiler are:
# daKD_gse_dutta, daOE_gse_dutta, dadaOE_gse_dutta, scOE_gse_dutta
daKD_gse_dutta@organism <- 'dm'
daOE_gse_dutta@organism <- 'dm'
dadaOE_gse_dutta@organism <- 'dm'
scOE_gse_dutta@organism <- 'dm'
daKD_gse_dutta_gtr   <- importCP(daKD_gse_dutta, type = 'gsea')
daOE_gse_dutta_gtr   <- importCP(daOE_gse_dutta, type = 'gsea')
dadaOE_gse_dutta_gtr <- importCP(dadaOE_gse_dutta, type = 'gsea')
scOE_gse_dutta_gtr   <- importCP(scOE_gse_dutta, type = 'gsea')

plotGSEA(scOE_gse_dutta_gtr,
         plot_type = "classic",
         show_pathway = scOE_gse_dutta_gtr$gsea_df$ID[1:2],
         #show_gene = scOE_gse_dutta_gtr$genelist[seq(1, by=10, length=2),'ID']
         show_gene = c('da', 'sc', 'emc', 'pros', 'poxn', 'Myo31DF', 'ck', 'nub', 'pdm2', 'esg', 'Dl'))

Layered heatmap to display NES across conditions and sets

Prepare data to plot heatmaps with NES and padj with the function below, gseCP_summarise. This is likely more complicated than it needs to be, as (now) all the GSEA objects have the same number of gene sets, but that is because I am running GSEA with pvalueCutoff = 1. If you run it with pvalueCutoff<1, non-significant gene sets do not even show up in the @result slot. In case this is happening, gseCP_summarise does a full join.

You run this separately for the NES or adjusted p values for “historical” reasons (I didn’t think it through, and had to restore to new syntax that now I don’t want to throw away :).

Now we can apply gseCP_summarise to the GSEA results from all conditions:

gse_list <- list(scOE_gse_dutta, dadaOE_gse_dutta, daOE_gse_dutta, daKD_gse_dutta)
conditions <- c('*scute*', '*da:da*', '*da*', '*da^RNAi^*')
# gse_list and conditions must have the same order
sets.as.factors <- c("ISC-only", "Progenitor", "EB-only", "EE-only",
                     "Absorptive", "EC-only", "Differentiation")
layerhm.df <- gseCP_summarise(dutta_gmx, gse_list, conditions, sets.as.factors,
                          cluster=TRUE, nsig.out = FALSE)

subt <- "for cell-type specific gene sets (Dutta et al., 2015)"
p <- layer.heatmap(layerhm.df, subt)
p + geom_hline(aes(yintercept=3.5), linewidth = 0.5) +
  theme(plot.margin = margin(r = 0))

2.2 RNAi screen phenotypes (Zeng et al., 2015)

To simplify later on the GSEA calls:

# common gene set
zeng_gmx <- import_from_gmx(paste0(getwd(),'/resources/Zeng.gmx'))
zeng_gmx$term <- str_replace_all(zeng_gmx$term, c(
  "Excess of.ISC.or.EB" = "Excess ISC/EB",
  "ISC to.EC"           = "ISC to EC",
  "ISC to.EE"           = "ISC to EE",
  "Overproliferation"   = "Over-proliferation" ) )
# common parameters
GSEAparams <- list(exponent = 1, minGSSize = 1, maxGSSize = 1000, eps = 0, pvalueCutoff = 1,
                   pAdjustMethod = "BH", TERM2GENE = zeng_gmx, TERM2NAME = NA, verbose = TRUE,
                   seed = FALSE, by = "fgsea")

Now we can simply define the rank file and apply GSEA for each condition:

# da knockdown
zrank <- make_degrank(DaKD_deg, mode='log2fc', key='ensemblGeneID')
daKD_gse_zeng <- do.call(
  GSEA, c(list(geneList=zrank), GSEAparams)
  )
# da overexpression
zrank <- make_degrank(DaOE_deg, mode='log2fc', key='ensemblGeneID')
daOE_gse_zeng <- do.call(
  GSEA, c(list(geneList=zrank), GSEAparams)
  )
# da:da overexpression
zrank <- make_degrank(DaDaOE_deg, mode='log2fc', key='ensemblGeneID')
dadaOE_gse_zeng <- do.call(
  GSEA, c(list(geneList=zrank), GSEAparams)
  )
# scute overexpression
zrank <- make_degrank(ScOE_deg, mode='log2fc', key='ensemblGeneID')
scOE_gse_zeng <- do.call(
  GSEA, c(list(geneList=zrank), GSEAparams)
  )

Layered heatmap

First to produce the dataframes for ggplot2:

gse_list <- list(scOE_gse_zeng, dadaOE_gse_zeng, daOE_gse_zeng, daKD_gse_zeng)
# `conditions` were defined further above
sets.as.factors <- c("Excess ISC/EB", "Over-proliferation", "ISC to EE", "ISC to EC", "Large nucleus", "ISC death")
layerhm.df <- gseCP_summarise(zeng_gmx, gse_list, conditions, sets.as.factors,
                              cluster=TRUE, nsig.out = FALSE)

subt <- "for RNAi-induced phenotypic class gene sets (Zeng et al., 2015)"
p <- layer.heatmap(layerhm.df, subt)
p + geom_hline(aes(yintercept=3.5), linewidth = 0.5) +
  theme(plot.margin = margin(r = 0))

GSEA plots

Prepare the data:

daKD_gse_zeng@organism <- 'dm'
daOE_gse_zeng@organism <- 'dm'
dadaOE_gse_zeng@organism <- 'dm'
scOE_gse_zeng@organism <- 'dm'
daKD_gse_zeng_gtr   <- importCP(daKD_gse_zeng, type = 'gsea')
daOE_gse_zeng_gtr   <- importCP(daOE_gse_zeng, type = 'gsea')
dadaOE_gse_zeng_gtr <- importCP(dadaOE_gse_zeng, type = 'gsea')
scOE_gse_zeng_gtr   <- importCP(scOE_gse_zeng, type = 'gsea')

Tentative plot:

plotGSEA(scOE_gse_zeng_gtr,
         plot_type = "classic",
         show_pathway = scOE_gse_zeng_gtr$gsea_df$ID[3],
         show_gene = c() )

2.3 Gene List Annotation for Drosophila (Hu et al., 2015)

Prepare the data for sub.group analyses

gladdir <- paste0(getwd(),'/resources/GLAD')
gladfiles <- paste(gladdir, list.files(gladdir), sep='/')
names(gladfiles) <- lapply (
  gladfiles,
  \(x) str_split(str_split(x, "GLAD_")[[1]][[2]], '_')[[1]][[1]] )
read.add <- function(csv, csvlist) {
  df <- read.csv(csv)
  df$term <- names(which(csvlist==csv))
  return(df)
}
datas <- lapply(gladfiles, \(x) read.add(x, gladfiles))
glad_dataset <- do.call(rbind , c(datas, make.row.names=FALSE))

Tidy up the strings for term names:

substitutions <- c(' signaling pathway$'='',
                   '\\.'=' ',
                   'Co '='Co-',
                   'Non '='Non-',
                   'RTK| RTK'='',
                   "DNA "="DNA-",
                   "TNF "="TNF-",
                   "TGF "="TGF-",
                   ' and '=' / ',
                   'transcription factor|Transcription factor'='TF',
                   ':'='/',
                   '^HLH'='bHLH',
                   '^Basic'='b',
                   '\\|'='')
glad_dataset <- glad_dataset %>%
  # remove 1-2-1 substitutions above
  mutate(across(c(term, Sub.group, Sub.sub.group),
                \(x) str_replace_all(x, substitutions))) %>%
  # remove trailing spaces
  mutate(across(c(term, Sub.group, Sub.sub.group),
                \(x) str_replace_all(x, '  $| $', ''))) %>%
  # start all strings with Uppercase
  mutate(across(c(term, Sub.group, Sub.sub.group),
                \(x) gsub(x, pattern="(^[[:lower:]])", replacement="\\U\\1", perl=TRUE))) %>%
  # substitute empty strings for NAs
  mutate_at(c('term', 'Sub.group', 'Sub.sub.group'), ~na_if(., ''))
View(glad_dataset %>% dplyr::select(term, Sub.group, Sub.sub.group) %>% distinct() %>% filter(term=='TF/DNA-binding'))

Establish the gene sets in tidy dataframes (subdivisions are in lists):

# major groups
glad_gmx <- glad_dataset %>%
  dplyr::select(term, FBgn) %>%
  rename(gene = FBgn)

# subgroups
glad_sub_gmx <- refine_glad_by(glad_dataset, 'Sub.group')
glad_sub2_gmx <- refine_glad_by(glad_dataset, 'Sub.sub.group')

GSEA calls, major lists:

# common gene set
#old_glad_gmx <- import_from_gmx(paste0(getwd(),'/resources/GLAD.gmx'))
# common parameters
GSEAparams <- list(
  exponent = 1,
  minGSSize = 1,
  maxGSSize = 5000,
  eps = 0,
  pvalueCutoff = 1,
  pAdjustMethod = "BH",
  TERM2NAME = NA,
  verbose = TRUE,
  seed = FALSE,
  by = "fgsea"
  )

Now we can simply define the rank file and apply GSEA for each condition:

# da knockdown
rank <- make_degrank(DaKD_deg, mode='log2fc', key='ensemblGeneID')
daKD_gse_glad <- do.call(
  GSEA, c(list(geneList=rank, nPermSimple = 1000, TERM2GENE = glad_gmx), GSEAparams)
  )

# da overexpression
rank <- make_degrank(DaOE_deg, mode='log2fc', key='ensemblGeneID')
daOE_gse_glad <- do.call(
  GSEA, c(list(geneList=rank, nPermSimple = 1000, TERM2GENE = glad_gmx), GSEAparams)
  )

# da:da overexpression
rank <- make_degrank(DaDaOE_deg, mode='log2fc', key='ensemblGeneID')
dadaOE_gse_glad <- do.call(
  GSEA, c(list(geneList=rank, nPermSimple = 1000, TERM2GENE = glad_gmx), GSEAparams)
  )

# scute overexpression
# rank <- make_degrank(ScOE_deg, mode='log2fc', key='ensemblGeneID')
# scOE_gse_glad <- do.call(
#   GSEA, c(list(geneList=rank, nPermSimple = 1000000, TERM2GENE = glad_gmx), GSEAparams)
#   )
# saveRDS(scOE_gse_glad, 'output/scOE_gse_glad.RDS')
scOE_gse_glad <- readRDS('output/scOE_gse_glad.RDS')

Layered heatmaps, main GLAD terms

First to produce the dataframes for ggplot2:

gse_list <- list(scOE_gse_glad, dadaOE_gse_glad, daOE_gse_glad, daKD_gse_glad)
# `conditions` were defined further above
sets.as.factors <- unique(glad_gmx$term)
layerhm.df <- gseCP_summarise(glad_gmx, gse_list, conditions, sets.as.factors,
                              cluster = TRUE, nsig.out = TRUE)

subt <- "for Gene List Annotation for Drosophila gene sets (Hu et al., 2015)"
p <- layer.heatmap(layerhm.df, subt)
p + geom_hline(aes(yintercept=3.5), linewidth = 0.5) +
  theme(plot.margin = margin(r = 100))

GSEA and layered heatmaps for GLAD sub.terms

Now let’s produce the gsea for the sub-terms of GLAD:

daKD_gse_gladsub   <- subglad_gsea(DaKD_deg, glad_sub_gmx)
daOE_gse_gladsub   <- subglad_gsea(DaOE_deg, glad_sub_gmx)
dadaOE_gse_gladsub <- subglad_gsea(DaDaOE_deg, glad_sub_gmx)
# scOE_gse_gladsub   <- subglad_gsea(ScOE_deg, glad_sub_gmx, perms=100000)
# saveRDS(scOE_gse_gladsub, 'output/scOE_gse_gladsub.RDS')
scOE_gse_gladsub <- readRDS('output/scOE_gse_gladsub.RDS')
scOE_gse_gladsub['Human Disease'] <- NULL

Now it is a matter of looping over the names of the GLAD terms with sub/sub-groups, and gather the NES/p.adj data for each condition and each term subdivision. For the 'Sub.group' subdivision:

gladsub_gsealists <- list(daKD_gse_gladsub,
                          daOE_gse_gladsub,
                          dadaOE_gse_gladsub,
                          scOE_gse_gladsub)

# 'loop' over the Sub.groups in `glad_sub_gmx`
# collect the corresponding GSEA results for each of the conditions
# gather the NES
lhm_gladsub_list <- lapply(
  1:length(glad_sub_gmx),
  \(x) gseCP_summarise(glad_sub_gmx[[x]],
                       # this is needed instead of the simpler `lapply(gladsub_gsealists, '[[', x)`
                       # because `gladsub_gsealists` and `glad_sub(2)_gmx` do not coincide in indexing anymore
                       # this can be corrected by running again the GSEA but scute data take a long time
                       # will do that overnight.
                       lapply(gladsub_gsealists, '[[', names(glad_sub_gmx)[[x]]),
                       conditions,
                       unique(glad_sub_gmx[[x]]$term),
                       cluster = TRUE,
                       nsig.out = TRUE)
  )

names(lhm_gladsub_list) <- names(glad_sub_gmx)

And plotting:

subt <- paste0("for Gene List Annotation for Drosophila gene subsets (Hu et al., 2015)",
               "\n(subsets of ", names(glad_sub_gmx), ")")
p_list <- lapply(
  1:length(glad_sub_gmx),
  \(x) if(nrow(lhm_gladsub_list[[x]])>0) {
    layer.heatmap(lhm_gladsub_list[[x]], subt[[x]]) +
    geom_hline(aes(yintercept=3.5), linewidth = 0.5) +
    theme(plot.margin = margin(r = 50))
    }
)

And to visualise:

p_list[[2]]

p_list[[1]]

## NULL

p_list[[3]]

p_list[[4]]

 p_list[[5]]

## NULL

p_list[[6]]

p_list[[7]]

GSEA and layered heatmaps for GLAD sub.sub.terms

Run GSEA:

daKD_gse_gladsub2   <- subglad_gsea(DaKD_deg, glad_sub2_gmx)
daOE_gse_gladsub2   <- subglad_gsea(DaOE_deg, glad_sub2_gmx)
dadaOE_gse_gladsub2 <- subglad_gsea(DaDaOE_deg, glad_sub2_gmx)
scOE_gse_gladsub2   <- subglad_gsea(ScOE_deg, glad_sub2_gmx)

gladsub2_gsealists <- list(daKD_gse_gladsub2,
                           daOE_gse_gladsub2,
                           dadaOE_gse_gladsub2,
                           scOE_gse_gladsub2)

lhm_gladsub2_list <- lapply(
  1:length(glad_sub2_gmx),
  \(x) gseCP_summarise(glad_sub2_gmx[[x]],
                       lapply(gladsub2_gsealists, '[[', names(glad_sub2_gmx)[[x]]),
                       conditions,
                       unique(glad_sub2_gmx[[x]]$term),
                       cluster = TRUE,
                       nsig.out = TRUE)
  )

names(lhm_gladsub2_list) <- names(glad_sub2_gmx)

And plotting:

subt <- paste0("for Gene List Annotation for Drosophila gene sub-subsets (Hu et al., 2015)",
               "\n(subsets of ", names(glad_sub_gmx), ")")
names(subt) <- names(glad_sub_gmx)
p2_list <- lapply(
  names(glad_sub2_gmx),
  \(x) if (nrow(lhm_gladsub2_list[[x]])>0) {
    layer.heatmap(lhm_gladsub2_list[[x]], subt[[x]]) +
    geom_hline(aes(yintercept=3.5), linewidth = 0.5) +
    theme(plot.margin = margin(r = 50))
    }
)

p2_list[[1]]

p2_list[[2]]

p2_list[[3]]

p2_list[[4]]

NAs are probably due to small groups where the members are simply not detected in the DEG set.

GSEA plots

For sub-sub-terms

for (x in 1:length(daKD_gse_gladsub2)) daKD_gse_gladsub2[[x]]@organism <- 'dm'
for (x in 1:length(daOE_gse_gladsub2)) daOE_gse_gladsub2[[x]]@organism <- 'dm'
for (x in 1:length(dadaOE_gse_gladsub2)) dadaOE_gse_gladsub2[[x]]@organism <- 'dm'
for (x in 1:length(scOE_gse_gladsub2)) scOE_gse_gladsub2[[x]]@organism <- 'dm'

daKD_gse_gladsub2_gtr <- list(
  importCP(daKD_gse_gladsub2$`Major signaling pathways`, type = 'gsea'),
  importCP(daKD_gse_gladsub2$Matrisome, type = 'gsea'),
  importCP(daKD_gse_gladsub2$`TF/DNA-binding`, type = 'gsea'),
  importCP(daKD_gse_gladsub2$Transporters, type = 'gsea')
  )
names(daKD_gse_gladsub2_gtr) <- names(daKD_gse_gladsub2)

daOE_gse_gladsub2_gtr <- list(
  importCP(daOE_gse_gladsub2$`Major signaling pathways`, type = 'gsea'),
  importCP(daOE_gse_gladsub2$Matrisome, type = 'gsea'),
  importCP(daOE_gse_gladsub2$`TF/DNA-binding`, type = 'gsea'),
  importCP(daOE_gse_gladsub2$Transporters, type = 'gsea')
  )
names(daOE_gse_gladsub2_gtr) <- names(daOE_gse_gladsub2)

dadaOE_gse_gladsub2_gtr <- list(
  importCP(dadaOE_gse_gladsub2$`Major signaling pathways`, type = 'gsea'),
  importCP(dadaOE_gse_gladsub2$Matrisome, type = 'gsea'),
  importCP(dadaOE_gse_gladsub2$`TF/DNA-binding`, type = 'gsea'),
  importCP(dadaOE_gse_gladsub2$Transporters, type = 'gsea')
  )
names(dadaOE_gse_gladsub2_gtr) <- names(dadaOE_gse_gladsub2)

scOE_gse_gladsub2_gtr <- list(
  importCP(scOE_gse_gladsub2$`Major signaling pathways`, type = 'gsea'),
  importCP(scOE_gse_gladsub2$Matrisome, type = 'gsea'),
  importCP(scOE_gse_gladsub2$`TF/DNA-binding`, type = 'gsea'),
  importCP(scOE_gse_gladsub2$Transporters, type = 'gsea')
  )
names(scOE_gse_gladsub2_gtr) <- names(scOE_gse_gladsub2)

plotGSEA(daKD_gse_gladsub2_gtr[[1]],
         plot_type = "classic",
         show_pathway = daKD_gse_gladsub2_gtr[[1]]$gsea_df$ID,
         show_gene = c('da', 'sc', 'emc', 'pros', 'esg', 'Dl'))

plotGSEA(scOE_gse_gladsub2_gtr[[1]],
         plot_type = "classic",
         show_pathway = scOE_gse_gladsub2_gtr[[1]]$gsea_df$ID,
         show_gene = c('da', 'sc', 'emc', 'pros', 'esg', 'Dl', 'N', 'H', 'Su(H)', 'E(spl)malpha-HLH'))

For sub-terms

for (x in 1:length(daKD_gse_gladsub)) daKD_gse_gladsub[[x]]@organism <- 'dm'
for (x in 1:length(daOE_gse_gladsub)) daOE_gse_gladsub[[x]]@organism <- 'dm'
for (x in 1:length(dadaOE_gse_gladsub)) dadaOE_gse_gladsub[[x]]@organism <- 'dm'
for (x in 1:length(scOE_gse_gladsub)) scOE_gse_gladsub[[x]]@organism <- 'dm'

daKD_gse_gladsub_gtr <- list(
  importCP(daKD_gse_gladsub$Kinases, type = 'gsea'),
  importCP(daKD_gse_gladsub$`Major signaling pathways`, type = 'gsea'),
  importCP(daKD_gse_gladsub$Matrisome, type = 'gsea'),
  importCP(daKD_gse_gladsub$Metabolic, type = 'gsea'),
  importCP(daKD_gse_gladsub$Phosphatases, type = 'gsea'),
  importCP(daKD_gse_gladsub$`TF/DNA-binding`, type = 'gsea'),
  importCP(daKD_gse_gladsub$Transporters, type = 'gsea')
  )
names(daKD_gse_gladsub_gtr) <- names(daKD_gse_gladsub)

daOE_gse_gladsub_gtr <- list(
  importCP(daOE_gse_gladsub$Kinases, type = 'gsea'),
  importCP(daOE_gse_gladsub$`Major signaling pathways`, type = 'gsea'),
  importCP(daOE_gse_gladsub$Matrisome, type = 'gsea'),
  importCP(daOE_gse_gladsub$Metabolic, type = 'gsea'),
  importCP(daOE_gse_gladsub$Phosphatases, type = 'gsea'),
  importCP(daOE_gse_gladsub$`TF/DNA-binding`, type = 'gsea'),
  importCP(daOE_gse_gladsub$Transporters, type = 'gsea')
  )
names(daOE_gse_gladsub_gtr) <- names(daOE_gse_gladsub)

dadaOE_gse_gladsub_gtr <- list(
  importCP(dadaOE_gse_gladsub$Kinases, type = 'gsea'),
  importCP(dadaOE_gse_gladsub$`Major signaling pathways`, type = 'gsea'),
  importCP(dadaOE_gse_gladsub$Matrisome, type = 'gsea'),
  importCP(dadaOE_gse_gladsub$Metabolic, type = 'gsea'),
  importCP(dadaOE_gse_gladsub$Phosphatases, type = 'gsea'),
  importCP(dadaOE_gse_gladsub$`TF/DNA-binding`, type = 'gsea'),
  importCP(dadaOE_gse_gladsub$Transporters, type = 'gsea')
  )
names(dadaOE_gse_gladsub_gtr) <- names(dadaOE_gse_gladsub)

scOE_gse_gladsub_gtr <- list(
  importCP(scOE_gse_gladsub$Kinases, type = 'gsea'),
  importCP(scOE_gse_gladsub$`Major signaling pathways`, type = 'gsea'),
  importCP(scOE_gse_gladsub$Matrisome, type = 'gsea'),
  importCP(scOE_gse_gladsub$Metabolic, type = 'gsea'),
  importCP(scOE_gse_gladsub$Phosphatases, type = 'gsea'),
  importCP(scOE_gse_gladsub$`TF/DNA-binding`, type = 'gsea'),
  importCP(scOE_gse_gladsub$Transporters, type = 'gsea')
  )
names(scOE_gse_gladsub_gtr) <- names(scOE_gse_gladsub)

plotGSEA(daKD_gse_gladsub_gtr[[1]],
         plot_type = "classic",
         show_pathway = daKD_gse_gladsub_gtr[[1]]$gsea_df$ID,
         show_gene = c('da', 'sc', 'emc', 'pros', 'esg', 'Dl'))

plotGSEA(scOE_gse_gladsub_gtr[[1]],
         plot_type = "classic",
         show_pathway = scOE_gse_gladsub_gtr[[1]]$gsea_df$ID,
         show_gene = c('da', 'sc', 'emc', 'pros', 'esg', 'Dl', 'N', 'H', 'Su(H)', 'E(spl)malpha-HLH'))

For main terms

daKD_gse_glad@organism <- 'dm'
daOE_gse_glad@organism <- 'dm'
dadaOE_gse_glad@organism <- 'dm'
scOE_gse_glad@organism <- 'dm'

daKD_gse_glad_gtr <- importCP(daKD_gse_glad, type = 'gsea')
daOE_gse_glad_gtr <- importCP(daOE_gse_glad, type = 'gsea')
dadaOE_gse_glad_gtr <- importCP(dadaOE_gse_glad, type = 'gsea')
scOE_gse_glad_gtr <- importCP(scOE_gse_glad, type = 'gsea')

plotGSEA(daKD_gse_glad_gtr,
         plot_type = "classic",
         show_pathway = daKD_gse_glad_gtr$gsea_df$ID[c(2,3,6)],
         show_gene = c('da', 'sc', 'emc', 'pros', 'esg', 'Dl'))

plotGSEA(scOE_gse_glad_gtr,
         plot_type = "classic",
         show_pathway = scOE_gse_glad_gtr$gsea_df$ID[c(2,3,6)],
         show_gene = c('da', 'sc', 'emc', 'pros', 'esg', 'Dl', 'N', 'H', 'Su(H)', 'E(spl)malpha-HLH'))